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Abstract: We propose an Adaptive Perfectly Matched Layer (APML) to 
be used in diffraction grating modeling. With a properly tailored co-ordinate 
stretching depending both on the incident field and on grating parameters, 
the APML may efficiently absorb diffracted orders near grazing angles (the 
so-called Wood’s anomalies). The new design is implemented in a finite 
element method (FEM) scheme and applied on a numerical example of a 
dielectric slit grating. Its performances are compared with classical PML 
with constant stretching coefficient. 
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1. Introduction 

Since their introduction by Berenger in H] for the time dependent Maxwell’s equations, Per¬ 
fectly Matched Layers (PMLs) have become a widely used technique in computational physics. 
The idea is to enclose the area of interest by surrounded layers which are absorbing and per¬ 
fectly reflectionless. Subsequent formulations (||2l O) allow to consider PMLs as a complex 
change of co-ordinates, which implies equivalent material properties (I?] |5l). A properly de¬ 
signed change of variables preserves the solution in the region of interest while introducing 
exponential decay at infinity. It is then obvious to truncate the problem to a finite domain, set 
a convenient boundary condition on the boundary of this domain and apply the FEM. In a 
scattering problem, the “traditional co-ordinate stretching” is frequency dependent, yielding a 
constant attenuation of propagating waves. However, these kinds of PMLs are inefficient for pe¬ 
riodic problems when dealing with grazing angles of diffracted orders, i.e. when the frequency 
is near a Wood’s anomaly (0121), leading to spurious reflexions and thus numerical pollution 
of the results. An important question in designing absorbing layers is thus the choice of their 
parameters : the PML thickness and the absorption coefficient. To this aim, adaptive formula¬ 
tion have been set up, most of them employing a posteriori error estimate (0|9l[TOl). In this 
paper, we propose Adaptive PMLs (APMLs) with a suitable co-ordinate stretching, depending 
both on incidence and grating parameters, capable of efficiently absorbing propagating waves 
with nearly grazing angles. In the first section, we detail our FEM formulation and expose the 
problems arising when dealing with Wood’s anomalies. Section 2 is dedicated to the mathemat¬ 
ical formulation used to determine PML parameters adapted to any diffraction orders. In the last 
section, we provide a numerical example of a dielectric slit grating showing the efficiency of 
our approach in comparison with classical PMLs. 

2. Setup of the problem and notations 

2.1. Governing equations and description of the structure 

We denote by x, y and z the unit vectors of an orthogonal co-ordinate system Oxyz. We deal 
with time-harmonic fields, so that the electric and magnetic fields are represented by complex 
vector fields E and H with a time-dependence in exp(— /tuf), which will be discarded in the 
sequel. Moreover, we denote ko = (o/c. 

The materials in this paper may be z-anisotropic, so the tensor fields of relative permittivity 
e and relative permeability jj. are of the following form : 
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where the coefficients Exx, ena, - fe are (possibly) complex valued functions of x and y, and 
where (resp. JIf) is the complex conjugate of £„ (resp. /r„). 

The grating is illuminated by an incident plane wave of wave vector defined by the angle 9o : 
k+ = a^x + P^y = A:+(sin Pqx — cos 9oy). Its electric (resp. magnetic) held is linearly polarized 






Fig. 1. Set up of the problem and notations 


along the z-axis, this is the so-called transverse electric or s-polarization case (resp. transverse 
magnetic or p-polarization case): 

E*’ = (resp. = A° e(''"^-'')z) (2) 

where A^ (resp. A® ) are arbitrary complex numbers and r = (x,y)^. 

The diffraction problem we are dealing with is to solve Maxwell’s equations in harmonic 
regime, i.e. to find the unique solution (E,H) of; 

fcurlE = iat/iofiH (3a) 

\curlH = —itueoeE, (3b) 

such that the diffracted field E‘* ;= E — E** and H'* := H — H** satisfies an Outgoing Wave 
Condition (OWC, ifTTI ) and where E and H are quasi-periodic with respect to the x co-ordinate. 

Under the aforementioned assumptions, the diffraction problem in non conical mounting 
can be separated in two fundamental scalar cases TE and TM. Thus we search for a z-linearly 
polarized electric (resp.magnetic) field E = e(x,y)z (resp. H = h{x,y)x). Denoting e and p. the 
2x2 matrices extracted from e and ji : 

e=( ^ ) and p = ( ) 
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the functions e and h are solution of similar differential equations ; 

•= div(|_-grad u)+klxu = 0, 

with 
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for the TE case, 
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for the TM case. 


The gratings we are dealing with are made of three regions (see Fig. [Til : 

• The superstrate (y > h^) which is supposed to be homogeneous, isotropic and lossless 

and characterized solely by its relative permittivity e+ and its relative permeability 
and we denote := where k{) := (ojc. 

• The substrate (y < 0) which is supposed to be homogeneous and isotropic and there¬ 
fore characterized by its relative permittivity e and its relative permeability and we 
denote k^ : = k^^Je^p^ 

• The groove region (0 < y < h^), embedded in the superstrate, which can be heteroge¬ 
neous, z-anisotropic and lossy, so that it is characterised by the tensor fields e® {xy) and 
p^{x,y). The periodicity of the grooves along Ox is denoted d. 


2.2. Implementation of the PMLs 

Transformation optics have recently unified various techniques in computational electromag¬ 
netics such as the treatment of open problems, helicoidal geometries or the design of invisibility 
cloaks (ID). These apparently different problems share the same concept of geometrical trans¬ 
formation, leading to equivalent material properties. A very simple and practical rule can be 
set up (jSl) : when changing the co-ordinate system, all you have to do is to replace the initial 
materials properties e and p by equivalent material properties and ps given by the following 
rule: ~ ~ 


£,=J-i£j-Tdet(J) and = J-Vj-^detd), (8) 

where J is the Jacobian matrix of the co-ordinate transformation consisting of the partial deriva¬ 
tives of the new co-ordinates with respect to the original ones is the transposed of its in¬ 
verse). 

In this framework, the most natural way to define PMLs is to consider them as maps on a 
complex space C^, which co-ordinate change leads to equivalent permittivity and permeability 
tensors. We detail here the different co-ordinates used in this paper. 

• {x,y,z) are the cartesian original co-ordinates. 

• {xs,ys,Zs) are the complex stretched co-ordinates. A suitable subspace F C is chosen 
(with three real dimensions) such that {xs,ys,Zs) are the complex valued co-ordinates of 
a point on F (e.g. x = 9le(xj), y = 5Ie(yi), z = iHe(Zi)). 

• {xc,yc,Zc) are three real co-ordinates corresponding to a real valued parametrization of 

Fed 

In this paper, we use rectangular PMLs ( ifTSl l absorbing in the y-direction and we choose a 
diagonal matrix J = diag (I, Sy (y), I), where Sy (y) is a complex-valued function defined by : 

ys{y) = f sr(y)dy'. (9) 

Jo 

The expression of the equivalent permittivity and permeability tensors are thus ; 
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Fig. 2. The basic cell used for the FEM computation of the diffracted field Mj- 


Note that the equivalent medium has the same impedance than the original one as e an /r are 
transformed in the same way, which guarantees that the PML is perfectly reflectionless. ~ 

We are now in position to define the so-called substituted field = (EijH^), solution of Eqs. 
Q with ^ and ^ = Xs- equals the field F in the region <y <y‘ (with y* = —h^ 
and y =W -F7!+), provided that iv(y) = 1 in this region (cf. Fig.|2l). The complex co-ordinate 
mapping ^(yc) is simply defined as the derivative of the stretching coefficient ij, (y) with respect 
to yc- With simple stretching functions, we can obtain a reliable criterion on the decaying of the 
fields. A classical choice is : 


v{y) = 


if y < y“ 
if y'’ < y < y' 
if y > y' 


( 11 ) 


where + i^"'^ are complex constants with > 0. 

In that case, the complex valued function y(yc) defined by Eq. (|9]l is explicitly given by : 


yiyc) = {yc 


/+r(Tc-/) 

if yc < y’’ 


if y* < yc < y 

y+c+(Fc-/) 

if yc > y' 


( 12 ) 


Let’s now study the behaviour of the field in the PMLs. In the substrate and the superstate, 
according to Bloch’s theorem, the diffracted field u‘^{x,y) can be written as : 


«"(x,y) = ^ «^(y)e-''- 
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Eventually, denoting 


= (15) 

and accounting for an OWC, we obtain the expressions of the different diffraction orders ; 

f M+(y) = for y > 

4{y) = { (16) 


“n (y) = ^ for y < 0. 


We focus on the bottom PML, as similar considerations can be conducted for the top PML. 
We consider a diffraction order propagating in the substrate, which can be expressed in the 
stretched co-ordinate system as : 


= u-{y{yc)) = 

The non oscillating part of this function is given by : 

t/,7(y) = t „exp + py K) , 

' _ n _ / _ n _ 

where j3,7 = j3„’ +/j3n' . For a propagating order we have j3„’ >0andj3„’ = 0, while for 

/ _ // _ / n 

an evanescent order j3„’ = 0 and )3„' > 0. It is thus sufficient to take ^ > 0 and -^ > 0 to 

ensure the exponential decay to zero of the field inside the PML if it was of infinite extent. But, 
of course, for practical purposes, the height of the PML is finite and have to be suitably chosen. 
For this, two pitfalls must be avoided : 

1. The PML is too small compared to the skin depth. As a consequence, the electromagnetic 
wave cannot be considered as vanishing : this limited PML is no longer reflectionless. 

2. The PML is much larger than the skin depth. In that case, a significant part of the PML 
is not useful, which gives rise to the resolution of linear systems of large dimensions. 

It then remains to derive the skin depth, , associated with the propagating order n. This 
characteristic depth is defined as the length for which the field is divided by e : 
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and we take the larger value among the /„ : 
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We set the height of the bottom PML region h =101. 



2.3. The FEM formulation 

Our FEM formulation is that described in mini. It relies on the fact that the diffraction 
problem can be rigorously treated as an equivalent radiation problem with sources inside the 
diffractive object. The expression of the source terms depends on the field ui solution of the 
annex problem of a simple interface, and is known in closed form. In this context, the role of the 
PMLs is to damp the field radiated by the sources. We set quasi-periodicity conditions on lateral 
boundaries and Neumann homogeneous conditions are imposed on the outward boundary of 
each PML. By so doing, the electromagnetic wave is not forced to vanish and the values of the 
computed field on these boundaries give valuable information on the absorbing efficiency of 
the PMLs. The cell defined in Fig. |2]is meshed using 2''‘* order Lagrange elements ( ifTSl l. In 
the numerical examples in the sequel, the maximum element size is set to Ao/(N,„A/iHe(er)), 
where Nm is an integer (between 6 and 10 is usually a good choice). The final algebraic system 
is solved using a direct solver (PARDISO). 

Note that we use this FEM formulation in this paper but the results of the sequel about PMLs 
are general and can be applied to any numerical method used for the modelling of diffraction 
gratings. 

2.4. Weakness of the classical PML for grazing diffracted angles 
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Fig. 3. Zero* transmitted order by a grating with a rectangular cross section (see parameters 
in text, t)art l2.4t for different values of : blue line, ^ = 1, correct damping; green 

line, ^ =0.1, underdamping; red line, ^ = 20, overdamping. 


Resuming part l2.2l we consider only the bottom PML, as analogous conclusions holds for the 
top PML. The efficiency of the classical PML fails for grazing diffracted angles, in other words 
when a given order appears/vanishes : this is the so-called Wood’s anomaly, well known in the 
theory of gratings. In mathematical words, there exists «o such that ~ 0. The skin depth 
of the PML then becomes very large. To compensate it, one should increase the value of 
but this gives rise to spurious numerical reflections due to an overdamping. For a fixed value of 
h^, if ^ is too weak, the absorption in the PMLs is insufficient and the wave reflects on the 
outward boundary of the PML. To illustrate these typical behaviours (cf. Fig.[3ll, we compute 
the diffracted held by a grating with a rectangular cross section of height h^ = 1.5pm and width 
















= 3[im with = 11.7, deposited on a substrate with permittivity = 2.25. The structure 
is illuminated by a p-polarized plane wave of wavelength Aq = 10pm and of angle of incidence 
00 = 10° in the air (e+ = 1). All materials are non magnetic = 1) and the periodicity of the 
grating is c/ = 4pm. We set = IQIq and ^ = 1. 

3. Construction of an adaptive PML 

To overcome the problems pointed out in the preceding section, we propose a co-ordinate 
stretching that rigorously treats the problem of Wood’s anomalies. The wavelengths “seen” 
by the system are very different depending on the order at stake : 

• if the diffracted angle 0„ is zero, the apparent wavelength A / cos 0„ is simply the incident 
wavelength, 

• if the diffracted angle is near ±;r/2 (grazing angle), the apparent wavelength A / cos 0„ is 
very large. 

Thus if a classical PML is adapted to one diffracted order, it will not be for another, and vice 
versa. The idea behind the APML is to deal with each and every order when progressing in the 
absorbing medium. 

Once again the development will be conducted only for the PML adapted to the substrate. 
We consider a real-valued co-ordinate mapping yd{y), the final complex-valued mapping is 
then yc(y) = C yd{y)^ with the complex constant with ’^ > 0 and ^ ’^ > 0, accounting 
for the damping of the PML medium. We start by transforming Eq. (fTSl l. n a function 

with integer argument into a function with real argument continuously interpolated between the 
imposed integer values. Indeed, the geometric transformations associated to the PML has to be 
continuous and differentiable in order to compute the Jacobian. For this purpose, we choose the 
parametrization: 


a{yd) = ao + ^^, (17) 

a Ao 

so that the application defined by fi^{yd)^ = — <x{yd)^ is continuous. Thus, the 

propagation constant of the transmitted order is given by (nAo). The key idea is to 

combine the complex stretching with a real non uniform contraction (given by the continuous 
function y{yd), Ea.(fT9ll). This contraction is chosen in such a way that for each order n there is a 
depth such that around this depth the apparent wavelength corresponding to the order in play 
is contracted to a value close to Aq. At that point of the PML, this order is perfectly absorbed 
thanks to the complex stretch. We thus eliminate first the orders with quasi normal diffracted 
angles at lowest depths up to grazing orders (near Wood’s anomalies) which are absorbed at 
greater depths. In mathematical words, the translation of previous considerations on the real 
contraction can be expressed as : 


exp[-/)3 (yrf)y(yrf)] =exp(-!koyrf) (18) 

The contraction y{yd) is thus given by : 
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The functiony(y^) has two poles, denoted^ — sin 0o). Wheny^ ^ = ±«Ao with 

n G N*, )3^(y^ j_) = j3^(±nAo) = = 0, i.e. we are on a Wood’s anomaly associated with the 





appearance/disappearance of the ±n* transmitted order. We now search for the nearest point to 
associated with a Wood’s anomaly, denoting : 


[</ D+=mm\y* -nXAo\ 

«*/ = min |y;5_ + «* Ao|, 

n^GN* ’ 

In a second step, we seek the point = n*Ao such that: 

n*/ D= min (D+,D^), (20) 

n*G{n* ,nf } 

To avoid the singular behaviour at y^/ = y^ j., we continue the graph of the function yf/(y) by 
a straight line tangent at y°, which equation is fo(yrf) = — y®) +y{y^d)^ where s{yd) = 

-^{yd) is the so-called stretching coefficient. The final change of co-ordinate is then given by : 

{ y{yd) foryrf <y° 

( 21 ) 

toiyd) foryrf >y0. 

Figure 0] shows an example of this co-ordinate mapping. 



Fig. 4. Example of a co-ordinate mapping y{yd) used for the APML (black solid line). The 
graph of yd{y) (blue solid line) is continued by a straight line tQ{yd) tangent at yJJ (red 
dashed line) to avoid the singular behaviour at yd = yj- 

Eventually, the complex stretch Sy used in Eq. (Soil is given by : 

^y{yd) = C^-^{yd)- ( 22 ) 

oyd 

Equipped with this mathematical formulation, we can tailor a layer that is doubly perfectly 
matched; 

• to a given medium, which is the aim of the PML technique, through Eq. (l8]l. 









to all diffraction orders, through the stretching coefficient ij,, which depends on the char¬ 
acteristics of the incident wave and on opto-geometric parameters of the grating. 


4. Numerical example 

We now apply the method described in the preceding parts to design an adapted bottom PML for 
the same example as in part 12.41 The parameters are the same, and we choose the wavelength 
of the incident plane wave close to the Wood’s anomaly related to the +1 transmitted order 
(Aq = 0.999y^ ^).Moreover, we set the length of the PML = l.ly^ ^ and choose absorption 
coefficients = 1+ i. For both cases (PML and APML), parameters are alike, the only 

difference being the complex stretch Sy. 
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(a) Classical PML (b) Adapted PML 

Fig. 5. Field maps of the logarithm of the norm of Ex and Ey for the dielectric slit grating 
at Ao = 0.999y^ y_ (same parameters as in part 12.4b . (a) : classical PML with inefficient 
damping of in the bottom PML. (b) : APML where the field is correctly damped in 
the bottom PML. For both cases the thickness of the PML is = 1. ly^ 

The field maps of the norm of H^, Ex and Ey are plotted in logarithmic scale on Fig. |5] for 
the case of a classical PML and our APML. We can observe that the field that is effectively 
computed is clearly damped in the bottom APML (leftmost on Fig. IHb)) whereas it is not in 
the standard case (leftmost on Fig. 13 a)), causing spurious reflections on the outer boundary. 
The fields Ex and Ey are deduced from thanks to Maxwell’s equations. The high values of 
Ey at the tip of the APML (rightmost on Fig. Hb)) are due to very high values of the optical 
equivalent properties of the APML medium (due to high values of Sy), which does not affect 
the accuracy of the computed field within the domain of interest. 

Another feature of our approach is that it efficiently absorbs the grazing diffraction order, as 
illustrated on Fig. |3; the -Fl transmitted order does not decrease in the standard PML (blue 
solid line), and reaches a high value at y = —h^, whereas the same order tends to zero as 
y —^ —h^ in the case of the adapted PML (blue dashed line). 

To further validate the accuracy of the method, we compare the diffraction efficiencies 































PML substrate 


Fig. 6. Modulus of the u„ for the three propagating orders with adapted (dashed lines) and 
classical PMLs (solid lines). Note that the classical PMLs are efficient for all orders except 
for the grazing one (n = 1) as expected. This drawback is bypassed when using the adaptive 
PML. 


computed by our FEM formulation with PML and APML to those obtained by another method. 
We choose the Rigorous Coupled Wave Analysis (RCWA), also known as the Fourier Modal 
Method (EMM, IfT^ ). For the chosen parameters, only the 0* order is propagative in reflexion 
and the orders —1,0 and +1 are non evanescent in transmission. We can also check the energy 
balance B = Rq + T_i +To + T^i since there is no lossy medium in our example. Results are 
reported in Table [T] and show a good agreement of the FEM with APML with the results 
from RCWA. On the contrary, if classical PML are used, the diffraction efficiencies are less 
accurate compared to those computed with RCWA. Checking the energy balance leads the 
same conclusions : the numerical result is perturbed by the reflection of the waves at the end 
of the PML if it is not adapted to the situation of nearly grazing diffracted orders. 



Ro 

T-i 

To 

T+i 

B 

RCWA 

0.1570 

0.3966 

0.1783 

0.2680 

0.9999 

FEM + APML 

0.1561 

0.3959 

0.1776 

0.2703 

0.9999 

FEM + PML 

0.1904 

0.4118 

0.1927 

0.2481 

1.0430 


Table 1. Diffraction efficiencies Rq, 71 1 , Tq and T+i of the four propagating orders, and 
energy balance B = Rq + T_i +Tq + T+i, computed by three methods : RCWA (line 1), 
FEM formulation with APML (line 2), FEM formulation with classical PML (line 3). 


Eventually, to illustrate the behaviour of the adaptive PML when the incident wavelength 
gets closer to a given Wood’s anomaly, we computed the mean value of the norm of along 
the outer boundary of the bottom PML 7 = {\H^{—h^)\)jt, when Aq = (1 + 10^")y^^ and 
Ao = (1 — 10 ^") 7 ^ for « = 1,2, ...10. The results are shown in Fig.|2] As the wavelength gets 
closer to 7 first increases but for 11 > 3, it decreases exponentially. However, in all cases, 
the value of 7 remains small enough to ensure the efficiency of the PMLs. 
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Fig. 7. Mean value of the norm of along the outer boundary of the bottom PML 7 = 
for Ao approaching the Wood’s anomaly ^ by inferior values (Aq = (1 — 
10 ^")yJ_l_, red squares) and by superior value (Aq = (1 + blue circles) as a 

function of n. 


5. Conclusion 

In this paper, we have proposed an adaptive PML that can treat rigorously Wood’s anomalies in 
numerical analysis of diffraction gratings. It is based on a complex-valued co-ordinate stretch¬ 
ing that deals with grazing diffracted orders, yielding an efficient absorption of the field inside 
the PML. We provided an example in the TM polarization case (but similar results hold for 
the TE case), illustrating the efficiency of our method. The value of the magnetic field on the 
outward boundary of the PML remains small enough to consider there is no spurious reflection. 
The formulation is used with the FEM but can be applied to others numerical methods. More¬ 
over, the generalization to the vectorial three-dimensional case is straightforward : the recipes 
given in part [3 do work irrespective of the dimension and whether the problem is vectorial. 
In addition, although the designed co-ordinate stretching is specific to the context of Wood’s 
anomalies, one can adapt this kind of non uniform PML to others wave problems by following 
the methodology exposed here. 




